{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Scalable GP Regression in 1D (w/ KISS-GP)\n",
    "\n",
    "## Introduction\n",
    "\n",
    "For 1D functions, SKI (or KISS-GP) is a great way to scale a GP up to very large datasets (100,000+ data points).\n",
    "Kernel interpolation for scalable structured Gaussian processes (KISS-GP) was introduced in this paper:\n",
    "http://proceedings.mlr.press/v37/wilson15.pdf\n",
    "\n",
    "SKI is asymptotically very fast (nearly linear), very precise (error decays cubically), and easy to use in GPyTorch!\n",
    "As you will see in this tutorial, it's really easy to apply SKI to an existing model. All you have to do is wrap your kernel module with a `GridInterpolationKernel`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import torch\n",
    "import gpytorch\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "# Make plots inline\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set up training data\n",
    "\n",
    "We'll learn a simple sinusoid, but with lots of training data points. At 1000 points, this is where scalable methods start to become useful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_x = torch.linspace(0, 1, 1000)\n",
    "train_y = torch.sin(train_x * (4 * math.pi) + torch.randn(train_x.size()) * 0.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set up the model\n",
    "\n",
    "The model should be somewhat similar to the `ExactGP` model in the [simple regression example](../01_Simple_GP_Regression/Simple_GP_Regression.ipynb).\n",
    "\n",
    "The only difference: we're wrapping our kernel module in a `GridInterpolationKernel`. This signals to GPyTorch that you want to approximate this kernel matrix with SKI.\n",
    "\n",
    "SKI has only one hyperparameter that you need to worry about: the grid size. For 1D functions, a good starting place is to use as many grid points as training points. (Don't worry - the grid points are really cheap to use!). You can use the `gpytorch.utils.grid.choose_grid_size` helper to get a good starting point.\n",
    "\n",
    "If you want, you can also explicitly determine the grid bounds of the SKI approximation using the `grid_bounds` argument. However, it's easier if you don't use this argument - then GPyTorch automatically chooses the best bounds for you."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "class GPRegressionModel(gpytorch.models.ExactGP):\n",
    "    def __init__(self, train_x, train_y, likelihood):\n",
    "        super(GPRegressionModel, self).__init__(train_x, train_y, likelihood)\n",
    "        \n",
    "        # SKI requires a grid size hyperparameter. This util can help with that. Here we are using a grid that has the same number of points as the training data (a ratio of 1.0). Performance can be sensitive to this parameter, so you may want to adjust it for your own problem on a validation set.\n",
    "        grid_size = gpytorch.utils.grid.choose_grid_size(train_x,1.0)\n",
    "        \n",
    "        self.mean_module = gpytorch.means.ConstantMean()\n",
    "        self.covar_module = gpytorch.kernels.GridInterpolationKernel(\n",
    "            gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel()),\n",
    "            grid_size=grid_size, num_dims=1,\n",
    "        )\n",
    "\n",
    "    def forward(self, x):\n",
    "        mean_x = self.mean_module(x)\n",
    "        covar_x = self.covar_module(x)\n",
    "        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)\n",
    "\n",
    "    \n",
    "likelihood = gpytorch.likelihoods.GaussianLikelihood()\n",
    "model = GPRegressionModel(train_x, train_y, likelihood)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Train the model hyperparameters\n",
    "\n",
    "Even with 1000 points, this model still trains fast! SKI scales (essentially) linearly with data - whereas standard GP inference scales quadratically (in GPyTorch.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter 1/30 - Loss: 1.142\n",
      "Iter 2/30 - Loss: 1.113\n",
      "Iter 3/30 - Loss: 1.085\n",
      "Iter 4/30 - Loss: 1.055\n",
      "Iter 5/30 - Loss: 1.024\n",
      "Iter 6/30 - Loss: 0.991\n",
      "Iter 7/30 - Loss: 0.958\n",
      "Iter 8/30 - Loss: 0.925\n",
      "Iter 9/30 - Loss: 0.888\n",
      "Iter 10/30 - Loss: 0.835\n",
      "Iter 11/30 - Loss: 0.753\n",
      "Iter 12/30 - Loss: 0.639\n",
      "Iter 13/30 - Loss: 0.513\n",
      "Iter 14/30 - Loss: 0.404\n",
      "Iter 15/30 - Loss: 0.320\n",
      "Iter 16/30 - Loss: 0.257\n",
      "Iter 17/30 - Loss: 0.202\n",
      "Iter 18/30 - Loss: 0.153\n",
      "Iter 19/30 - Loss: 0.106\n",
      "Iter 20/30 - Loss: 0.061\n",
      "Iter 21/30 - Loss: 0.015\n",
      "Iter 22/30 - Loss: -0.036\n",
      "Iter 23/30 - Loss: -0.077\n",
      "Iter 24/30 - Loss: -0.123\n",
      "Iter 25/30 - Loss: -0.164\n",
      "Iter 26/30 - Loss: -0.203\n",
      "Iter 27/30 - Loss: -0.245\n",
      "Iter 28/30 - Loss: -0.281\n",
      "Iter 29/30 - Loss: -0.310\n",
      "Iter 30/30 - Loss: -0.352\n"
     ]
    }
   ],
   "source": [
    "# Find optimal model hyperparameters\n",
    "model.train()\n",
    "likelihood.train()\n",
    "\n",
    "# Use the adam optimizer\n",
    "optimizer = torch.optim.Adam([\n",
    "    {'params': model.parameters()},  # Includes GaussianLikelihood parameters\n",
    "], lr=0.1)\n",
    "\n",
    "# \"Loss\" for GPs - the marginal log likelihood\n",
    "mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)\n",
    "\n",
    "training_iterations = 30\n",
    "for i in range(training_iterations):\n",
    "    optimizer.zero_grad()\n",
    "    output = model(train_x)\n",
    "    loss = -mll(output, train_y)\n",
    "    loss.backward()\n",
    "    print('Iter %d/%d - Loss: %.3f' % (i + 1, training_iterations, loss.item()))\n",
    "    optimizer.step()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Make predictions\n",
    "\n",
    "SKI is especially well-suited for predictions. It can comnpute predictive means in constant time, and with LOVE enabled (see [this notebook](../05_Scalable_GP_Regression_Multidimensional/KISSGP_Deep_Kernel_Regression_With_LOVE_Fast_Variances_CUDA.ipynb)), predictive variances are also constant time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQAAAADDCAYAAABtec/IAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvFvnyVgAAIABJREFUeJztnXt8k+W9wL9JW5qCNKEV2sqtLReFKVCIKHMybIs6mDAqCkeZbmPAmWy6c/CC4vGyqTgdmzh3q0XdxpjVnipyZLDRbiKCQqAVhVagrYJYaKGk5ZZekpw/nry5vuktaZI2z/fz4UOa95Ined/n9/6e31Vjt9uRSCTRiTbcA5BIJOFDCgCJJIqRAkAiiWKkAJBIohgpACSSKCY20BMYjcZcx8uZJpPpoUDPJ5FIQkdAGoBj8t9mMpm2AZONRuPk4AxLIpGEAk2w4gCMRmOVyWQaFZSTSSSSkBAUG4DRaHwQWBaMc0kkktARTA3gTWCJyWQyq21fuXKlDDmUSMLEs88+q1F7PyAjoLLmN5lM+4BqYCnwnL/9n3zyyQ7PWVdXx5AhQwIZVo8T6WOM9PFB5I8x0scHnR/j448/7ndboEuAXCDJ8dqAEAISiaSXEKgAyAcyjUbjUgCTyVQU+JAkEkmoCGgJ4Fjv5wdpLJIooa2tjbNnz3L27FkiNRvVZrPR1NQU7mG0i/cYNRoN8fHxpKamEhvbuakdcCCQRNJVTpw4gV6vJzk5GY1G1TYVdlpbW4mLiwv3MNrFe4x2ux2z2cyJEycYNmxYp84hQ4ElIae5uZnExMSInfy9FY1Gg8FgoLm5udPHSAEgCTl2uz0iJn9ZWRllZWU9/jlms5ni4uIe/xwQQqAryyopACQRTW1tLbm5uZw4caLb5ygrK6OgoICSkhIKCgqorhbOKr1eT1FRz9utDQaD6ueUlZUxbtw4iouLKS4uZs2aNc6xqdHetu4ibQCSiGb16tXs3LmTZ555hhdffLHLx5vNZp5//nk2bNjgfO+OO+5gw4YNJCUltXNkcBk0aJDPe1lZWWRkZJCXl+d8b9asWWzevNln3+rqatatW8fTTz8d1HFJASCJSAwGAxaLxfl3fn4++fn56HQ6zGbVYFNVioqKyM7O9nhv0KBBlJSUMGXKFMrKyigpKaG8vJzFixezd+9eAPbu3cv8+fMpLS0lKSmJjIwMampqKCoqIiMjg8svv5wtW7awYcMGli9fzooVKwA89s/IyGDdunVMmjSJffv2dfp7K0/60tJSALKzsykvL6empoaysjL0ej2lpaVYrVZmzpxJZmZmp38Pb+QSQBKRVFRUsGDBAhISEgBISEhg4cKFVFZWdvlcjY2NfrdlZWWRk5PDpEmTWLduHeXl5ZSWlnLDDTfw6KOPMmXKFOfkz87OZtCgQTz99NPcddddznPk5eWRmZnps/+qVauYN28eOTk5ZGRkdGnMmZmZJCUlkZSUxFtvvUV2djYZGRlkZWX5bAsEKQAkEUlaWhqJiYk0Nzej0+mcnoPU1NQunSc7O9v5VFeoqakhJyfH4z1lOTBv3jwWL17MmjVraGlpQa/Xk5WV5dQiDAaDx7nXrFnDlClTnO95799VzGYzmZmZrFmzBr1ez6RJk5zvg1gKKNsmTpzosa07yCWAJGKpq6tjyZIlLF68mHXr1nXLEJiZmckDDzxAQUEBGRkZlJeX89JLLzm3m81mjyWAorLfcMMNzJw5k3Xr1jmfvooKbjabMRgMzJ8/n1WrVjmFwlNPPeWx/4oVK3jrrbeYNGmS89isrCznZ5eVlVFTU+P0ENTU1DjHpnxeY2Mj1dXVnDlzBrPZTE1NjXNbQ0MD1dXV1NTUeJy3KwQtG7AjVq5caZfJQKEh0sd35MgRRo4cGdGBNr0xEEjhyJEjjB492vn3448/7jcbUC4BJJIoRgoAiSSKkQJAIolipACQSKIYKQAkkihGCgCJJIqRAkAiiWKkAJD0acrKypg2bZpH2m91dbXPe9FKMFqDLXW8HCVbg0m6ik4XH5TzWCzqRTCysrKckYC//e1vAZEboMTVRzvBaA22zWQyKcVBczs6RiIJNXq93u+26upqCgoKKC4upqyszPn3K6+8QnV1NSUlJcyaNYuSkhJWrVoVwlGHhkCXAJmI0uAgSoJ3Py9REpVYLM1B+dcReXl5FBQU+MTje2fweWfa5eTkYDAYyMnJCSjpJlIJtCqwe0XgyUBhYMORSHqGnJwc7rjjDo/MPQW9Xk9mZiYZGRmsWbOGSZMmMXz4cI4ePYrZbFYt5tFXCEo2oKND0D5HhyC/1NXVdXiu3iBlI32MkT4+m82G1WoNyWeVl5fz8ssvM3z4cCZNmsSwYcPYs2cPZWVl7NmzhyeeeIL8/Hyys7NJT09nxIgRHDlyhFOnTnHkyBE2bdpEdXU1hw4dorq6mj179jhTdMONv9/QZrN1aq5BkLIBjUbjgyaTyW9LMJDZgKEk0scnswGDQ0RkAxqNxqXK5JdGQImkdxEML8AvjEZjldFoPBOkMUkkkhARqBFwG9B3LSQSSR9HRgJKJFGMFAASSRQjBYBE0g1C2e6rJ5FVgSVh5zf/qgro+J/cMKrd7WVlZezdu9dZm7+0tLTTHXaUwKDy8nJn8w9wtfty7+rTG5ECQNKnUWsN1tknt9lspqGhgZycHNU2Yn0hQlAKAEmfpqioyCf8d8WKFVRXV1NaWkpGRgaNjY3o9XrWrFnDihUrKC0t5YknnmDv3r3U1NRQUlLCo48+yq5duzCbzT7tvpRzKS3BGhoaPM719NNPO3v7KWOZNGmSxzHhykyUNgBJ1KG08frhD39ITk4ORUVFqkk/SpJQTk4OkydPBlBt9+WdUKR2rjVr1rB48WLy8vLIzs72OSZcSA1A0qeZP38+99xzj8d7JSUlAM4OP2azOeCkH/eEIlBfHijLCKWTkPcx4UAKAEmfxmAweLQGa2xsZNKkSTz11FMUFRWRlJREXl4eNTU11NTUOFttlZeX09TU5GwFtm/fPsrKylTbfXm3BPM+l3Lc888/73zqex/j3nMwlMjWYN0g0scY6eOTyUDBISKSgSQSSe9FCgCJJIqJKgFgs9k539wW7mFEPRqNhlAtPaMNu92ORqOq7avSp42ANpudQ3XnqG20UH+2hdPnW2iz2UjUxTEiKYERSf0ZPiiBfrFRJQfDTnx8PE1NTSQnJ3fpZpW0j91ux2w2Ex/f+UrLfVYA1J1tprSynvpzvgUjmyytfPpVK59+1YQuLoacKwaTeemAMIwyOklNTeWLL76gsbExYjUBm82GVhvZDwbvMWo0GuLj40lNTe30OfqcAGi12viw5gwff9m5m8vSauXdT04wYZie6zKTiI2J7IveF4iNjWXgwIER7amIdE8KBGeMfUoAnLW08Xb5V5gvtnb52P1fNlJrtnDT14YwqH+/HhidRBJ59JnH3TlLG291c/Ir1J9rprislqYAziGR9CaCIgAcZcHDxrnmNorLv6IxCBP3Qksbmz45QXNraMpWSyThJBi9AXOBPwLtJ2X3EOea23irzDX5mxrqWb/6fhY9/EsA1j32I+qOVdPW4jIGamP7YWtrxTDkMsx1x4mJiycmJobbfvozil58gnue/xMD+sUwZ0IaWq20UgeT2tpaFixYwOnTp6muriYhIYFHHnmERx99lLi4OFpbW9HpdIwcOZKvvvqKkpISJkyYEO5h91kC1gAchUGrgzCWLtPSZmNjea2H2r9tw++p+dTEL3/0HZ76bg5fVVV4TH64FFtbDGDHXHccAGtrMy2WC/z12ftpvnCOl/77Tj40lTH1GzM4ceJEaL9UH8Nqs3OyycLHH39MUlISo0ePZvfu3VRVVWG327lw4QKPPvooIEJbASwWC5999hlnz57luuuuIzk5mWnTprF//35yc3PlNQkivdYGYLfb2VZZR8OFFgAenpPFA9/6GrveFd3JLjaZsdsUNX4E8DBwEKgHLgAW4CugDPhvwOUGbG228MKPb+XTfR+RNeVqeeN1kdraWr55Qzav/t97jJ98Dd/8xnVcc801XLhwwU83mwGAuu+6tbWV8+fPU1ZWxjXXXMOOHTvIzMykf//+lJaW9uj3iAZC6gUIZmuw/V+dZ/+XZwE4e+YUlw7L4ETNZ157XQ/8DJjh9t45xNfWAWmOf5MQAuLXwEtAk3PvM6frmTp1KgBjxoxh06ZNXHbZZZ0aY7gIZ2swq83O4vse4qNdOzH9xy1YW1tU9roMuAG4zvHvSsAOHAEOAJ8C/wA+8DhKcevabDYAZs2axdGjR3vke0R6ezUIzhhDKgA667PsaL9jDReoNJ9jwIABNDXU84cHFnHe3OC11wPAaiAGuAhsBP6CuLHaEAJgEKKn6SPA14GnHcfdA/zN53NbW1u5+eab0Wg01NTUdCngItSEw4etNxhotlicf/tO/njgUeAhwD2LrQWhjF7u+JcHPAa869j3gN/PHDFiBCCWDcEm0uMAIPAx9rolwFlLG1sP1mG323l4ThY/v3OG1+TXA28BzyEm/7NAKvAfwGbE5AexBKhF3GTXAdnAvwADsAEhFNSx2+2kp6cH82v1eppbrfzPH95sZ4/pwMcIARCHuBYrgW8AiYhlwATgDmANQgub7TimAKE1qJOdnS2XaN0kGL0B54v/jPODMJ52sdrs/P3Tk1xstfLwnCzafJ4wXwNMwHcAMzAHodo30TH/QgiB+wAbQhvIpz0lSafTha2QQyRxocXKW+W1vPvmepWt/YDfA+8hnu4HEQJ3NvALhJrfjNACPkFoXvcDo4HfIK7FYsSy4HrVzy8tLXXaBiRdIxhegCKTyTTIZDIVBWNA7fFhTQMnzwpVr63N2+efCZQgbpx9CNV+k8ceGo2GlJGjSRioJy5exyWDkh3vu/8MLwK3IpYNSxznGOgzFp1Ox8KFC6msrAz8i/VibDY7Q5IH8b3rRzsNsC7igDeA/0RM8CeALGCnx14xMTHMnTuXhIQEt3frgXuB8cAWxHLtn8Dt7YzFJoVyF+k1S4BjDRcoO9rotPbjEec/GHGTpCBukusAUWhRo9US209HXLyOr03L4f4/bGTS9Jtpa2kmrp8OjVaL3W7z+rS3EUaqeuBmoBhvTcBisfDvf/+7B75p76G2tpbJ115Hv4T+KltjEU/zucBpYBrwJNDizABMSEhg4cKFVFVVUVhYyI033siyZcu8qtwcQWgLLyFsCIXACtpD3dMgUSPiBUBtbS3ZObkUfXAQO3aVBJ8BwP8BYxBP/jzE+l5gt9nAbuOZt/dSuWe701Vot9s5c/K42I4QFJ58hLhpTwK5wK98xnbixAmeeeaZ4HzRXsioUaOo3L+PC03e1ugY4K8ITeoMMBNxbQR2u534+Hiam5tJTEx0GlMLCwtZu3Ythw4dYtQo97gyG/AThIEW4JfAC37HJe0znSfiBcDq1avZtfMD3nntNxyvqsTqofrHIlTMqYhYpFkIN58nitB4+LV/kDVjNnHxOgDi4nVk3TCbKTlznILAkypgHmKN+hNgmc8e+fn5Uad2GgwGdDqd0x3niRb4E0JVbwRuRMRaeGKz2ViyZAknT5702ZaWlubnKf5LYCHietyHP0Pt4cOH0el06HTiOtfW1kojoR8iVgAoN1l+fj52u51d7xbywo9v9drrN4hJr6jqvjdT8mUjeORP/wQgMWkwCQMucaj/8bS1NJORNph+djVftcIuYKnb533T9zOSkykuLo6am6z9NOv/Ae5EGF5vQhhlXcTExLBw4UJ27drF2rVrKSz0thsIJk6cyLJly9i8ebPXE70QIVwUQ+1CvyPZvXs3IB4iO3fujGptzR8RKwAqKiqYm3eb82ntywKEcckCfBs47LE1JlasI21Wq5j4cTFcP/pSLo2xsHTpUj7Y8T5Lly7Fev4M/9r8NjU1Ndx62+1+ikD8GeFWjAP+F2FwdHH69Gnuu+++qLnJXvjrO8TGqaVMZyP89zbEUuwjj60ajQar1UpiYmKH/mtlOZCdnc3MmTO9tr6DiN4EeA1h8/Fl6tSpzoeIzWaLSm2tIyK2HoA+eQgnzrbQ2mxBo41xC+sFMQHzHa9/Cuz2OPbSoel895E1fLj5Tc6bTzE1fRBZww30i9XyxhuuJ87atWudr9PS0kgy6P2otSDcieMRwqYQYR9w1Rc8fFgIoPz8fOeN1huiybrKiUYLf371VRUXbApi3a9FGPtKPLb279+fqVOncvnll3dZS6qrqyM9PZ0pU6awa9cuamtrsdvXIjw+P0YYba9FLNk80Wg06HQ6Ll68iEajYc6cOR7XPdqJKAFQW1tLXl4ecXFx3P3Ybzi0X6iPg1LSaKj90rFXHPA6InjkTUQioienjn/Or5ffSly/eGrrTnOJrnNfs66ujmXLlvHZZ59x/Phxqqur3daiNoRqux8wIoJYnlI9z8KFC3n22Wc7+a17DwaDwU/EnRYx+VMR8RQ/89njwoULfPjhh2zZsgXoXFi4grJM8P38nwIZCC/BZsR1Oetx7OjRo6mqqiImJgar1cqhQ4ciOoIz1ETUEmD16tWUlZWxe/duln/7Gs6cFNl6rskPIrLvaoSbb4nqeeLidXz9prkcOHiw05MfXGrnli1b+OSTT1SKKzYBP3C8fgyYqHoed8t2X2LD1l0eRlQXq4AchA3mDoSw9GTu3LkBx0xUVFQwZ84ct2WaFWED+BgYC/zW55jDhw9js9mcgryiokIuA9yICAHgbvBrn9mItV8r4sI3+uwR6zDujR8+hBHDhgY0rpkzZ7Jo0SKvd0sRxsA4hG3Ady2cn5+PXq/vU0bB0+daqLkYT/n2v9Pa7P4Uvg54HDHpFwHq3zclJSVgoZiWlkZKSorXMu0cwih4HvguQkvzJDU11SnMZQCXJxEhABTJHhMT085eQ4BXHa8fwXvdDzAldy6//ss7LFmyhPr6zquY/igsLKR/fxHkcskll7htWYkwOk5A3Py+tLa29hmjoN1up/iDA7x0/11e7tJ4RJy+knOxzbUlYQApqWnMnz+f9PR0VXdfd6irq1MpJX4I4RYEEXbsaaQ9ceIEzc2iJoTFYumzGlp3iAgbgCLZ24/gehkR8bcNkSziyZScOTz0zAvcclUqmrm+rrqu4r3ePHfOPb7gAvA94H1EttpGvAWS8pTqC0bBj79sYkP+C9R8agKNFpyRk6uAK4AKhOHPhc3axhefB7/tdWFhIbW1taxcuZJ33nmHixcvOrasQ7gdb0Mkc30DdyOtO/n5+fz5z3/utdcjmESEBgAuS+9V1+Uy7poZ4kZz8gNEYo8Z+D4id9wTW8tFbho3JGiNJioqKliwYIFXfLo7OxGBKTGIm09dlirhrr1V5dTr9Vw7JsUV5++c/FchNCGAHyJi/V20tjQ7A3GCTVpaGomJiTQ3i89wXfOlwFHgGrwFkoISh7B9+/Y+tUTrLhEjAAoLC6msrKRi93YqPvq3242WgSvscznwpc+x46+ezuaN/0t8XHtLiK7hfZNptVqVGIHHELHqVyISVzxRC3ftbTzzylsq72oRqn8cwvC202ePkSNHOgNxeoK6ujqWLFnC9u3bWbp0qWOpZkbYAKwI4fR1n+OsViuvv/4606ZN69VLNJstOA1VIkYAKHhGmSlhpQMRIb8bfPbXaLS89fbbPVLL3/0mW7JkiUpAihIiDCLTzTNnvbm5GY1GE7T1b6g52nCBog1/UtnyE0T49TFEfIQvX3zxBdOnT++xsSkem+nTp/PHP/6RCxcuOLbsQKQZK0JKvdSY1Wr1CA7S6/U9Ntae4Ej9+aCcJ+IEwD3PveaWmLMCkQNeC/xIdX+73cbEUT1Toku5ySZMmMDatWvZuHGjihawBZEtOBDvhKHY2Fg2bdrkN9w1ktEbDIy9LEklxXckIgQXxDU5iz9CsexRlmqeBuSfIewS4xChyYKYmBg0Gg2jRo1yLu2U42677bYeH2uwqK2tZd6sXObOnRvwEiaiBIDBYOA3/32nw9I8BVegzWLAu+SXIFTr6/YTYH6KcEMtQGQOCtra2rj33nt75Vrzb1t3+dmyFpGBWYiopqTOokWLQrLsGTduHIWFhV4G5GbEPWNDGGlFvIbVamXevHl8+eWXTuOhctz69et7TXzA/zz5c6oOlFNWVhbwEiZiBICn1f0SRLRfP4TP/e9+jwvV+to3CMWdY8DPHa9/i3tswJEjR9ixY0evSlHV6/XkfTNLZctsRH5/E0LoqTN+/HjOnvWvGQSTiooKhg4dquJC3oW4d2KBV1CMtMXFxU67zowZM5zH9QZjrfIQWv/aK873As1viBgB4Ln2/y0izrscVw64JzfOnsuyZctCtr5WD0IRiFJUv0aonWMRJa186S1PmLXF28maMRtw96joENWSAM0TeAf8DLhEVE0aP348Y8aMCdmyJy0tjVmzZmG3251eh379FAG8ChExOhnvIiJKQRer1YpOp4tYY617KnNFRQUjRoz02K5UU+qu4IoYAVBZWcnIkSMR0WR3IVRqJffbRWLSEFKGDichlnbTSXsCxVWZmprqrE0PSpvmNoSXAkSgkmcUouJ+qqysjOj89EMnz3HcfJEjH3+Ep7v1QUSAzX6wv+hz3Plz4ol/8OBBNm7cGFJB526sXbZsGcnJyQwcOBBxDymp3E8gisb40tzczJ133hmRxlollXnUqFFkZGRw9OgXHtutVmtAUZbBaA02H+F/mWwymZ7r7nnS0tKwWIYDv3O8cy/gqvPfLz6B5b9az+Xjr+KOq4fRLzb0sksRNnq93mPyfv75545X/0IkKN2GSB92haUq7qe3336bu+66y+mCevFF38kULppbrbx/5DRPfTfbK+IvA5e1fznCzSaqKGVedTUXG2o5+dVxrFYrCQkJzJ07N6TJUO4PgT/96U9eCUPbEEuAHwB/QOQsuIiJicFms5GQkMDLL78cgtF2Dr1e74xe9IdOp2Pw4MEBCa6AZpHSFNTRHswcSJNQnW4gJ0/+GmFNfx1x0Vy0NF/k18tv5b6brwrL5Hdnx44dHkFCOp3OWZ9eqP8XEUkx33CuMZX/LRZLxOanD05OYnnO5SrVkdYilgB/QbjZBHabjepP9jDrphudKni4VWnFK+DJA4iiMdmIfAEXVqsVu90ecddC8Uoo941vUNX3sVgmcOzYMbZu3drtzwl0Ji1APP1B1OTKbWdfv4gf3YBI8qlBrfQWQPaseXz2mXf3n9CTkpLiESTU0tLCgAFKa7GjCD80wItYrUKNtlqtJCcnM2fOHKfg0Gg0QcmSCwYnGi2sfPUfKrURb3H8a0TNHnP9jGzq6+s94iXCqUorAVyeNOAqIPIrIFn1WK1WG/Zr4TT0rRcl1hUvhadWMxJhJ9sFjA6oKUqgSwADnv459V/Wgb8c8Pfff5+nnnqKjRunA8PwruOv0WrBbictaQBarbZLueQ9gdls5tixY9x5553ceeedzJ49m4qKCrc9nkOonFkId5RQLU+fPs0777yDRqNBq9Vis9k4ePBg0L9TV2PcbXY77xxoQBvvHfacgNPwx2OolVzbXrqN+Ph4Z0GUVatWAR3n+/dUHP6YMWP8qM7rEfkbOYjrs9hnj7y8POrr61mwYAGrV6/ukfF1hDIXtm7d2s7E/iXi2mwAjrBly5Zu3z8R0RpsyJAhjm1twOce2+ITBjBs7JVcO+lKzplPRUy7prffftv5OjY21ssPfRGxFHgDETTzJi5FSXg8FK/H4cOHGTFiRNCThbryO+07aqb+TCOv/exe4nX9sdmstFguIqzo6QhvjG+uPYin5meffdat69IT17KyspKVK1eyceNGlQn0I0RBlx8gUrnf89ja2trKyy+/zJ49e3jppZeoq6tj/fr1IV3OKHOhpaXFWcQkMTGRixYLrS0tiD6X8xEGzgcBuPnmm4HutUcLVACYgSTHawOiAHy3qKurc35hd5ovnudYRTkHPnrPz5Hhp7KykhkzZvDFF+4W2jeBfyMu2BP485uHw2jmzpkLLXxUc4ZtG37Psc/2A6CNiUF08VFU/h+hGP68ueOOOyLKdaYsAVpa1Aq9HgaeQUQK/hGRzu3ab+PGjc7XxcXFgCh9fv58cMJuO8u6des83M1NTYpGHINLI1sNiII5I0eO5M0322vL5p9AbQCFuJKvM3FPCO/qiQoLOXLkCBOvv9mjbPfU3Fs44KFeRx5paWl+shDvQ0yc5YjsOV/CaTSz2eyMSB3Mf900ziPk12a14gpoehn4UPX4UAb8dAXFLagetPULRLzG5bTX/1FBiRMIpXGwqqpKJbwZhG3sKoSdzJUS379/fyZMmNCtzwpIAJhMpn0ARqMxFzArf3eXtLQ0dP0H0NbS7KzsM3roYIYPjex23CDKWKenpzN8+HC3d/cjJlIs3ip0YmIimzdvDovRTIlD+Oe+w6x8dStXTsvxMv4tRKyVT+FK+RXE9otn2g03kp6eHtKAn66g5HAoE8kzpbsFl5FZKfTqn3BECKqHNyfhijZdgXvzm0CWjsHoDZhvMpm2mUymjup5dUhtbS2f7NzG5Jw5fP+J3zLo0sGcPnE80NOGBCWdefLkyW4eAXAZz67H3QXV1NTEuHHjOHDgQMir1K5evZoPdu7kF88+Q2LSYC4ZlOzm+kvEldT0EN45GNNuymPrpreprKyMyMnvjrIccBUNUXgfsQRQNByX9jZ8+HAGDRoECPtGODQ0JbzZU6t8EiEEShDdr12cOXOm258VMZGAIG7Mi+ea6Bev47MPSzCfPuWIDuw9FBYWMnjwYLd33N1nzyPal8OUKVOYNm0aO3bs4Nprrw1JVKBHsxWbjZ3/9zoPfOtrfLj5Dbe9ngHSEDn+rzrfHTf1m0ybvZDYlsawx2F0BcW25MtDwFeImgGuTNNjx45x5swZxo4dy4cffhhSDU3RzC6//HKOHz/uFh6fhcsOc5/HMYFqJxFxJdW6AG1/528RGSzTHsr3cEUGKvwF8dRJQVHj9u7d65z0J06cID09vce/oxIkE68TKnFcvI6kFPeQ5W8i7BWtiKYrrlDgit3vsevd19n1XmmPjjHYKLYlz16DIATzjx2vVyPczy4OHTrE1KlTee2113pU03EPC1fCfocOdb9sQCpQAAAP8ElEQVQmMQgtRTEAHvA4PlDtJCIEgHf5La1W26uytBS8v4fnk2c5ws15DzBJ9XiLxdKjhSnS0tLQxvenxVHVt7XZQsNJZYnVH1HaDEQa9icex2q0Wm6cdUuvuA7e+O81+BailkMi/tycPV0nYPXq1c5sUSVCtKbGvZbivYjU+C9wr20AMH/+/IC1k4gQAO7lt2L79XPWcY+E0NKu4F1GzPOm+wSRnhqDqFyr/tP35A13ztLGwepjXDtrAZ6ZfiBU/1EIn79vEIzdZiN9WFqvuA5qTJyo3sNBaAGNiJqT/+GztafqBHSuFP5IXIa/exC+f8Gi7/2AX/3qVwFrJxEhAMDlullfWEx6ejrp6ekREVraVdwz0xYtWsTQoUPdrNCPI3y31+Kvx31P3XBWm52/vfcJ+z8oYde7r+OZ6XcdosxXK6LoaqvHsbFx/Rg6dFivug7eFBYW8tFHH7nlbCjU4roWv8N7KRAfH8/gwYPZvn17UMfTvrao8DtcxVc2e2xZ/9orjBmjnt3YFSKiLDi4Mrrq6uo81Mze1sfNXSIXFBQwYMAAN03gLCIEdQtCsr8LHARwhgb3VGDQ9sOn2JD/Ana7nfj+l9B8QSlznoBIvNIiohbLfY5ta23h+PEvOX36VFDHFGomTpzo5aFRWIfId5iLqEGZiyIgm5ubqa+v58UXX6SgoCCo43nvvfewWCzodDqVKL7bEZ2vzXgb/kAY/+6/X73uRFeIGA2gr5Kbm8vo0aPdglK2Iow68YibTchgJfKrJ5Y8iXo935o43Bns45r8IOLixyKWKOq9DhUsFkuvMMa2h9lsZvx4Nd//EoS7Nhu1qM1ga2arV692GoF9J/+luCL+HsA9B0Or1aLRaDrVYbkzSAHQw7z33nscOXLEq5LQCoRRR2ky6sJmswVV1d59sJrU9Mu58us5xPbzrpC7ELEGbkYkyqiFzwrcC5r0Zmpqati3b59K38d6XAlCq/EXuelZuarreK/9fc+nQeQppCBCydd5bBXFZ7RBu0ekAOhh1GvWnUWstcG9yahWq2X9+vWcPn06KHEBpopqZud+k2OHPqH+y89pa3HPkhuPKJsN4onXfhCn1WrlzTff7LVGQG8qKytVXIPvIgKE4hHZg74lxRUDb3dR1v7+z3E/8C1EFOYi1JrgWK3WgGoAuCMFQA+jVrNO8C/gJURzjb8AA7DZbCxatIgdO3YEXO1Vr9fzjazxnG2ox263c/JoldvWS4D/RRiY/oKolOOf+Ph4Ro8eTW5ut8o9RCT+XYMrcPV99HUNxsfHs3nzZp/3u/K5iYmJHmq/6+FwLcIbA0IjU4+CDWbdAikAQkBdXR02m01lrfcQouzZVYiJ6HLNKQFQCQkJXdYGRo/2lxOvUIDo6fcpIuDHPzqdjtbWVrKzsz2y5foCEydOZNGiRaSlpbm9q5R3v4BYEvyXxzHNzc0eqeBdxWAw+JQeE4LIgKiEFYtI9PFfch0ImiYmBUAIUFxQnhFeIG6yW4AzwDzUjHB2u71L2kBF7Vn+63cbyZoxG41WzbX0IOIGbwJudYzBE602hqFDh7Jo0aJe6YrtLIWFhRQUFKgI2DLgbsfr5xHWeBeBRKf6XwK8gvD778ZftyUQkX++Haq6jxQAIWLixIkqpapAqJu3IaIEH0Gtv71ywymW3/379/vs09JmY+vBk2yrrCPBkET59r9jt3mruPfiKlf2fURbbV9sNis33HADBQUFzq5IkZ74E3yKEPaZGOBvqGUNbt++vcsVntWWAKKk/DyEy28B7nEY3vEBSUlJQdXEpAAIEQaDwatsmDsluJqLrkOsBT3RaDSMHDmSpqYm7r77bo9tJ5ssvG76kkMnz9HUUM+rjy0j88qrufQy96CX5YjiniBKZRerjiR5cAogCp9GC9XV1SoGQRCxGq8jQoU3IdrTuygoKHDG7/vT0j7++GOGDBlCaWkpubm56PV6ryXAUwgjbAti8n/ucbzdbmfZsmXs3r2bZcuWtRPR2D00gbo1OsvKlSvtTz6p3rLZnbq6uogp++WP7oyxtraW++67j3fffdeP8QmEUXA5wgL8HeCDDs/71t4vqKo/j91up6mhnjU/+g4Xmrzzw/8TEX4MIqusfaOfO8EuVaYQadf5iiuuUEniAlEN+T1EM9TDwM2I+rfqeP9eWVlZVFRUYDAYaGpq4jvf+Q6bNm2itbUVofE9jdD+5gOeT3atVkt1dbXf9X5nf8PHH3+cZ599Vq1ijdQAQoXSWcj/5AfxJHgHEQhSgncJa3cGDU5l6Ojx7KuswW638/CcLH5+5wyvya9BBJIok385/ib/hOtvYm7ebc7Q1N6UhBUMJk6cSF5enlcqN4jCG7cg3KRjEJV4jT7He1d41ul06HQ6p9ZnNpux2WwUFxc7Jv99iMlvQ1xnX7XeZrNxxRVXBOkbqiMFQAhROgvdeuutzJ8/39lZyEUbkIeIAotHBIT8HN/EHbBcvMBXVRVs2/B7Hp6TRVurdxDPSKAUEekH4ob7Hf745IN/MiTZ4PRz96YkrGBQWFjIhg0buPTSS1W21iFSpbcCQxABOp6GQbvdztatW6mrqyM3N5fNmzer5B2A6HvxMvCC4+8fIpYZnmg0GoYNG9bjAlgKgBCiVA3661//yvr165k5cybjxo3z2ksp+qCkDz+KeDp4RqZdPNfkrJ3gO/m/hyhHNgMRRnoLrtBSTzTaGJKGpFFdVeWRyNRXLf8dYTabGThwoJdrEOAc8G3gNUT8xEaEQTXJuYfFYmHq1Kns2LGDWbNmcfToUa9zzEBclx8iNIsf4l50xR273c6sWbN6XAAHJRnIaDRODrQeYDRSWFioogUo/A44gigtrjTn2IpwS5Wo7D8KsY68HdEME0Swz38ibAp+sNuYe8ts0tLSPCz9vS0JK1goufjq16UN4T05hsjNfxDx+65BWPL9FUidgMg1UAqQmBBuxoOqe8+YMYNhw0KTfRmwBuAoCNq9msQSqqqqmDNnjp+t/0A8+dciAlRuQhRePoZYk/4b8SQqQwiLZxGTvwERRjqf9ia/NjaWu36wmIZT9cH4Kn2Kjz76qJ2Q3ccQdoAtCA/Bk4hKvZsQhtwHEIIiH3GtPkZM/lZESvg0/E1+gLFjx1JQUBAS12vAGoDJZNpmNBr9m0Ul7aIYB/1zDGEcfBLxtLkXkbM+zGu/RoQBsQghONpvEtFPl8AHez7mqjFq61SJ0v6tpaXFmartyV5EzP71CGPe9YglghpfIYTFb1BLt3YnPT09pEuvkNYD6Ez7op5qGRVMgj3GY8eOMXz4cE6dOqVSwVbhDCJL7XmEgU+PeProEdrBdtrL5vMmOTmZFL0ubG3WIv06u7d/q6qq4sCBAzQ2Nqrs+T4wHfgaMBpxbUYClyHW+5sRGkD7pKWlsWnTJqdbL1RzJSJag3V3v3ASzDEqseULFiwgNTWVxYsXc/fdd/sJHGoDqlTe909MbBxxOh0aNGi1MXx77jwsTafD/juH+/M7wj3mf8GCBfzjH/+gpaXFjyv3AN4FO7vCLbfcwpVXXtnl4wL9DTsUAEajcanK29WOluCSIOK+5hs7dizTp09n8eLF3H777Rw9elRFDfVPctpwWpubiY2L47JR47j7f9aS1L8fN38theRL+vXE8Ps0yrVRhPSrr75Ka2trl+oDxMbGYrfbSUlJ4cyZM6SkpPDGG2+wbt26kJSFVx1TRzsEo+GHpOu4CwN3X3BLm41tlXVU1fvvV3f+/Hmf0lfjUgcyfcylvaqmfySiXJdgeknC6XEJhhdgvvjPOD8I45F0QL9YLbOuTOW2yUMZavBu5+3L4EvimTluCLnjhsjJL/EhGF6AIoTpWRJCUvU68rIu44vTF6g4cZaLrVYsrTaaW63E2WLJSh/E2JRLGNRfqvsS/0RMVWBJ9xiZ3J+RyZ5BKyJJJMnPERKJC6kTSiRRjBQAEkkUIwWARBLFSAEgkUQxUgBIJFGMFAASSRQjBYBEEsVIASCRRDFSAEgkUYwUABJJFCMFgEQSxUgBIJFEMVIASCRRjBQAEkkUIwWARBLFSAEgkUQxUgBIJFGMFAASSRQTcEkwt7Lho0wm00OBnk8ikYSOgDQAR1/AbY7S4ZmOvyUSSS8h0CVAJqBM+mrH3xKJpJcQ0BLAq2nIZKDddqayN2BoiPTxQeSPMdLHBxHUG9BoNE4G9plMpn3t7Sd7A4aOSB8fRP4YI318EDm9AXOlAVAi6X0E3BvQaDQuNZlMzzle58qmoRJJ7yEYXoBfGI3GKqPReCZIY5JIJCEiUCPgNmBQkMYikUhCjIwElEiiGCkAJJIoRgoAiSSKkQJAIolipACQSKIYKQAkkihGCgCJJIqRAkAiiWKkAJBIohgpACSSKEYKAIkkipECQCKJYqQAkEiiGCkAJJIoRgoAiSSKkQJAIolipACQSKIYKQAkkigmGK3BlMYgM2VlYImkdxGMoqC3OWoDTnb0B5BIJL2EYBQFVcqAZ3bUGEQikUQWweoM9CCwrKP9Hn/88WB8nEQiCRIau90elBMZjcY3gSUmkynym6pJJBIgwNZgyprfofpXA0uB54I7RIlE0lME2hosF1DW/QZgTzAGJZFIQkNASwCj0WgAbnf8OcVkMnVoB5BIJJFD0GwAkvBgNBrnA2ZgstKk1c9+D7a3XRL5GI3Gyf48bZ29D7wJihegu3Q06O5+qRCOT7GPjApHEJSbDWab0WjM9HeDOOI1ZhIG+0wnfsPJQCaAyWQqCvHwlDF09j7M7Khbdk/huIZ/BEapbOvUfaBG2EKB3QcNmL2DiDraHgHjywW2OW6ITLeIyFCyAHFjgjDChmMMfunkNXzYMfEzwxFI1sn7sNqxvTpcwW7K5/vZ3O37IJy5AB0NOtw3d0efn+n2XrXj71BjABrc/k723sHxNNjm/X6IaPc3dDxZ9wCYTKbnwhRI1pn77BeO/yM12K3D+8Af4RQAHQ26218qSLT7+SaTKd9NHZwMmEI1sC6SFMbP7ugaXg0kG43GyY5gsnDQ0XXeh3jyn/Har08gswEDxKES7gvTk8GMa4IbgNPuG8P89O8sp5XfzqERRBQOT5cZWA28bDQaw6HpdUS790F7hFMAdDTobn+pINHZz88NYxZkIa6lRyaOvAzHTQtiXT3fYaxMCsP6taPf8DSuda0ZoRGEmo7GuBRY7TAOLgEiRki5XWfV+6AzhFMAdHTzdvtLBYmOxofRaFyqWI3DYQR0e3LmAmY3LaTEsb3IzbJuUDlFT9PRb1jktj1cgWQdXmcFx28ZllB3h3Zk9NKSlOvs7z7okLDGATieTNW4uVeMRuNek8k0xd/2SBmf48d+E7EuTMKVFi1xo5PXuAG4OlyaVCfG+KBje1K43IA9hQwEkkiiGGkElEiiGCkAJJIoRgoAiSSKkQJAIolipACQSKIYKQAkkihGCgCJJIr5fzxiyRw2Z+ZKAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 288x216 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Put model & likelihood into eval mode\n",
    "model.eval()\n",
    "likelihood.eval()\n",
    "\n",
    "# Initalize plot\n",
    "f, ax = plt.subplots(1, 1, figsize=(4, 3))\n",
    "\n",
    "# The gpytorch.settings.fast_pred_var flag activates LOVE (for fast variances)\n",
    "# See https://arxiv.org/abs/1803.06058\n",
    "with torch.no_grad(), gpytorch.settings.fast_pred_var():\n",
    "    test_x = torch.linspace(0, 1, 51)\n",
    "    prediction = likelihood(model(test_x))\n",
    "    mean = prediction.mean\n",
    "    # Get lower and upper predictive bounds\n",
    "    lower, upper = prediction.confidence_region()\n",
    "\n",
    "# Plot the training data as black stars\n",
    "ax.plot(train_x.detach().numpy(), train_y.detach().numpy(), 'k*')\n",
    "# Plot predictive means as blue line\n",
    "ax.plot(test_x.detach().numpy(), mean.detach().numpy(), 'b')\n",
    "# Plot confidence bounds as lightly shaded region\n",
    "ax.fill_between(test_x.detach().numpy(), lower.detach().numpy(), upper.detach().numpy(), alpha=0.5)\n",
    "ax.set_ylim([-3, 3])\n",
    "ax.legend(['Observed Data', 'Mean', 'Confidence'])\n",
    "\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
